Relevance of angular momentum conservation in mesoscale hydrodynamics simulations 
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The angular momentum is conserved in fluids with a few exceptions such as ferrofluids. However it 
can be violated locally in fluid simulations to reduce computational costs. The effects of this violation 
are investigated using a particle-based simulation method, multi-particle collision dynamics, which 
can switch on or off angular-momentum conservation. To this end, we study circular Couette 
flows between concentric and eccentric cylinders, where non-physical torques due to the lack of the 
angular-momentum conservation are found whereas the velocity field is not affected. In addition, 
in simulations of fluids with different viscosities in contact and star polymers in solvent, incorrect 
angular velocities occur. These results quantitatively agree with the theoretical predictions based 
on the macroscopic stress tensor. 

PACS numbers: 02.70.-c,47.11.-j,66.20.+d 



I. INTRODUCTION 

In simulations of the hydrodynamic behavior of com- 
plex fluids, one is faced with the challenge of bridging 
the gap between the mesoscopic length and time scales 
of the solute and the atomic scales of the solvent. As 
these length scales typically differ by orders of magni- 
tude, a full treatment on a microscopic level is prohibited 
by the huge number of involved particles and the large 
necessary time range. Moreover, one is often only inter- 
ested in the dynamics of the colloidal particles, while the 
microscopic details of the solvent that mediates the hy- 
drodynamic interactions are rather unimportant. Thus, 
a coarse-grained mesoscopic fluid model is required that 
is sufficiently simple to be tractable but still captures the 
correct hydrodynamic behavior. 

Various mesoscopic approaches have been proposed in 
the last decades. A large number of physical solvent 
molecules is represented by one model fluid particle at 
a time, reducing the number of degrees of freedoms con- 
siderably. Lattice methods, such as lattice gas automata 
(LGA) 1] and lattice-Boltzmann methods (LB) @, ||, 
generally suffer from the lack of Galilean invariance. 
Moreover, it is difficult to incorporate complex and de- 
formable boundaries that play important roles in the 
phase separation of two fluids [1, [j] and the dynamics 
of vesicles and cells [5] . In particle-based techniques such 
as dissipative particle dynamics (DPD) [1,0, [Hi or multi- 
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particle positions and velocities are continuous variables 
that are updated at discrete times. Coupling to so- 
lute particles as well as moving boundaries can be easily 
treated. MPC needs less computational time compared 
to other particle based methods such as DPD, thus al- 
lowing simulations of larger systems. 

In this article, we will focus on MPC, which has been 
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applied to various systems such as colloids [l|| ES H] , 
polymers [|, O S HH, Hip, [H , membranes jH, EH , 
ternary amphophilic fluids [2jj , and chemical reaction sys- 
tems 28]. The MPC method naturally contains ther- 
mal fluctuations. Hybrid simulations combining a MPC 
fluid with molecular dynamics (MD) of solute particles 
are easily possible. The algorithm is constructed in such 
way that mass, energy and translational momentum are 
locally conserved, which is essential for correct hydro- 
dynamic behavior. However, the angular momentum is 
not conserved in the most widespread version of MPC, 
which is often called stochastic-rotation dynamics (SRD). 
Here we denote it as MPC-SR. The consequences of this 
fact have not yet been investigated and are the subject 
of this paper. In order to clarify the effects of angular- 
momentum conservation, we mainly use the Andersen- 
thermostat version of MPC, denoted MPC-AT, where 
angular momentum conserving and non-conserving algo- 
rithms are available [29] . We also checked that the same 
quantitative dependence appears in the original MPC-SR 
method. The main conclusion is, that simulations that 
do not conserve angular momentum can lead to quan- 
titatively and even qualitatively incorrect results, when 
the boundary conditions on walls are given by forces, flu- 
ids with different viscosities are in contact, or finite-sized 
objects rotate in fluids. 

The rest of this paper is organized as follows: In Sec. 
IIT1 we briefly discuss the effect of the non-conservation 
of angular momentum on the stress tensor. Note, that 
while in the MPC fluid the non-conservation of angu- 
lar momentum is an artifact of the simulation method, 
there are also real fluids, where angular momentum is not 
conserved. For example, in ferrofluids asymmetric stress 
arises [3(1 [3l| when the rotation of the suspended parti- 
cles is impeded by external helds. In Sec. IIII1 the algo- 
rithms for the angular-momentum conserving and non- 
conserving versions of MPC-AT are described. A sim- 
ple geometry to study rotating fluids is the flow between 
rotating coaxial cylinders, also called circular Couette 
flow. The simulation results for the angular-momentum 
conserving and non-conserving methods are compared in 
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Sec. lIVI In particular, binary fluid and branched polymer 
systems are investigated in Sec. II V CI and D, respectively. 
In Sec. |Vj we lift the restriction of coaxiality and study 
rotating eccentric cylinders. This geometry is of practical 
importance in journal bearings and microfluidic devices 
using rotating colloids Finally, we summarize our 

results in Sec. IVH 



II. MACROSCOPIC FLUID DYNAMICS 

In conventional viscous fluids that do conserve angular 
momentum, the viscous stress tensor has to be symmet- 
ric, i. e. a a p = op a . This symmetry is required by the 
fact that there is no stress expected in a uniformly ro- 
tating fluid (rigid body rotation) [jjjj, or alternatively, 
by the conservation of angular momentum [34[ ■ On the 
other hand, for a fluid without conservation of angular 
momentum, the above argument is no longer valid and 
we have to consider in general an asymmetric tensor. 

Here, we consider Newtonian fluids, i. e. the stress is 
proportional to the strain rate, so that the a a p are linear 
functions of the derivatives dv a /dxp [33[. We decompose 
the stress tensor in its symmetric and asymmetric parts. 
Then, the viscous stress is given by 



tap = A(V • v)5 aP 

_ ( dva_ dv/3 \ . / dva_ _ dvp 
\dxf3 dx a J \dx0 dx a 



(1) 



where a,/3(z {x,y,z}. Here, A is the second viscosity 
coefficient, and fj and fj are the symmetric and asymmet- 
ric components of the viscosity, respectively. The last 
term in Eq. |T]) is linear in the vorticity V x v, and does 
not conserve angular momentum. Thus, the last term 
vanishes [i. e. fj — 0) in angular- momentum-conserving 
systems. 

The equation of velocity evolution is given by 

= -VP + (A + fj - rj)V(V • v) + (fj + rj)V 2 v, (2) 

where D/Dt is Lagrange's derivative and P is the pres- 
sure. When a fluid is incompressible, this is the normal 
Navier-Stokes equation with viscosity r\ = fj + fj. This 
is consistent with the usual definition of the shear vis- 
cosity r\ = o'xy/'j in simple shear flow with the velocity 
field v = jye x . Since the equations of continuity and 
velocity evolution are of the same form, the negligence 
of angular-momentum conservation does not modify the 
velocity field of fluids when the boundary conditions are 
given by velocities. However, it generates an additional 
torque, so that the velocity field can be changed when 
the boundary condition is given by forces. In cylindrical 
coordinates (r, 9, z) , the azimuthal stress is given by 



. .rd(vglr) „ vg 

<r r e = (fj + fj) V +23—- 
or r 



(3) 



angular velocity = vg/r. The second term is the addi- 
tional stress from the negligence of angular-momentum 
conservation and is proportional to il. 

When a fluid is compressible and the fluid density is 
not constant, the bulk viscosity is not negligible. The 
bulk viscosity without angular momentum conservation 
is given by A + 2rj/3 instead of A + 2ry/3. In angular- 
momentum-conserving fluids, these two values coincide 
because -q = fj. Thus, the effects of the angular- 
momentum conservation are not negligible when the 
torque on objects or the bulk viscosity is significant in 
fluid systems. Eqs. are general and can be applied 

to MPC methods and other model fluids, which do not 
conserve angular momentum. We explain the effects of 
the torque quantitatively using MPC- AT in the following 
sections. 



III. SIMULATION METHOD 

A. Multi-Particle Collision Dynamics 

MPC is one of the particle-based methods to simulate 
hydrodynamic behavior accompanied by thermal fluctu- 
ations. A fluid is described by point-like particles of mass 
m. The MPC algorithm consists of alternating streaming 
and collision steps. In the streaming step, the particles 
move ballistically, r,;(i + At) — i*j(i) + ViAt, where At 
is the time interval between collisions. Subsequently, the 
particles are sorted into the cells of a cubic lattice with 
lattice constant a that is randomly shifted before each 
collision step to ensure Galilean invariance [l(| ■ The col- 
lision step then mimics the simultaneous interaction of 
all particles within each cell by assigning the particles 
new velocities. There are several versions of the collision 
procedures and each version can switch on or off angular- 
momentum conservation [2^, [35| . We call the versions of 
methods with or without angular-momentum conserva- 
tion '+a' or '—a', respectively. In the original version 
(MPC-SR), the relative particle velocities with respect 
to the mean velocity in a cell are rotated by a fixed an- 
gle (p around an axis, which is chosen randomly for each 
cell In MPC- AT— a, the velocities of the particles are 
updated by [H HI] 



j £ cell 



ran 
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(4) 



where A^ c is the number of particles in a cell, and ve- 
locities vjj an are chosen from a Maxwell-Boltzmann dis- 
tribution. The center-of-mass velocity of each cell is 
conserved, and the temperature is constant in MPC- AT. 
In MPC-AT+a, the velocities of the particles are updated 

by m 
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The first term is the stress of the angular-momentum- 
conserving fluid, which depends on the derivative of the 
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where II is the moment-of-inertia tensor of the particles 
in the cell. The relative position is c = r,j — where 
rj? is the center-of-mass of all particles in the cell. 



B. Boundary Conditions 

In order to simulate no-slip boundary conditions, the 
following technique has been developed for —a fluids in 
Ref. In the streaming step, the fluid particles are 

scattered with a bounce-back rule on surfaces. In the 
collision step, in collision cells crossing a boundary with 
N c < n — (N c ), a virtual particle with mass m{n — N c ) 
and velocity v wa ii + v ran / (n — N c ) is inserted to calculate 
vj?, where v wa n is the velocity of the boundary wall. This 
algorithm keeps the slip on a boundary small [ill ]. 

We have tested some algorithms for +a methods, 
where the position of the virtual particle is now impor- 
tant. One possibility (denoted 'cen') is to locate it at the 
center of the cell. For a simple geometry like a cylinder, 
more sophisticated ways are available, e. g. by putting a 
virtual particle slightly inside boundary walls, which can 
reduce slip. A more direct way to estimate the velocities 
inside a wall is to distribute explicit particles inside the 
wall. Watari et al. [24] proposed a boundary algorithm, 
where particles freely enter inside objects and velocities of 
inside particles are updated to v ran . However, this allows 
flows to penetrate through a small object, when there is 
a pressure difference around the object. To prevent flow 
penetration, we employ the bounce-back rule. Particles 
are randomly distributed inside the cylinder wall with 
depth \J~2a from the surface with the same density as the 
outside fluid. Before collision steps, the velocity is up- 
dated to v wa rj + v ran . The position of the wall particles 
are updated by renewal of the random uniform distribu- 
tion [36|. In the Couette flow simulations, the velocity 
field is theoretically known. Thus, v wa u is extrapolated 
for a wall-particle position in most of the simulations (de- 
noted 'w-gra'). This explicit-particle boundary algorithm 
can be applied to other particle-based methods such as 
DPD. We employ 'w-gra' and 'cen' algorithms for coaxial 
systems (Sec. IIV|) and eccentric cylinders (Sec. IV)) . re- 
spectively We show the comparison of these two bound- 
ary algorithms and 'w-con' algorithm for Couette flow 
in Sec. IIVB1 In 'w-con', explicit wall particles with 
the constant angular velocity f^bd are employed so that 
Vwaii = fibd^ee. 



C. Viscosity 

The shear viscosity is calculated from u X y/i = V = 
fj + fj in simple shear flow with v = jyex . The viscosity 
of MPC consists of two contributions, r\ = 77km + rj co i; 
the kinetic viscosity 77km and the collision viscosity r) co \ 
result from the momentum transfer due to particle dis- 
placements and collisions, respectively. The viscosity of 
MPC- AT— a with large mean number density n, is given 
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FIG. 1: (Color online) Dependence of the viscosity 77 of MPC- 
AT-a (•, x) and MPC-AT+a (o, □) on the time step At for 
n = 10 in two dimensions. The inset shows the viscosity 
difference A7? co i = Jfcoi — jfcol of MPC- AT— a. 
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where d and k B Q are the spatial dimension, and thermal 
energy, respectively. The viscosity of MPC-AT+a with 
large n can be calculated similarly, and is found to be 



Vkm 
Vcol 



nk B QAt 



n- (d + 2)/4 



m(n — 7/5) 
24a d - 2 At ' 



(8) 
(9) 



The derivation and the correction terms for small n for 
+a versions of the MPC family will be reported else- 
where [35] ]. In two-dimensional systems, the angular- 
momentum constraint does not change the kinetic vis- 
cosity for large n, i. e. — ^kin- ^ e a ^ so calculate 
the viscosity from simulations for simple shear flow with 
Lees-Edwards boundary conditions [37]. Fig. [1] shows 
that the theoretical and numerical results are in very 
good agreement. 

The symmetric and asymmetric components of shear 
viscosity fj and 77 are calculated from the shear stress 
Vyx/i — fj — fj- Since the kinetic stress is symmetric in x 
and y, i. e. a^™ = cr^ n , the kinetic viscosity has no asym- 
metric component 77ki n = 0. The collision procedure of 
MPC- AT— a does not conserve the angular momentum. 
The molecular chaos assumption gives cr™ 1 = 0, because 
(v y (x)) = before and after the collisions. Thus, the 
viscosities are 



Vcd / 2- 



(10) 



This viscosity relation holds for all —a versions of MPC 
and DPD in Refs. [29, The numerical simulation of 
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FIG. 2: (Color online) Torque T in a rotating MPC-AT-a 
fluid with uniform angular velocity Qo for n = 10 and At = 
0.1. Symbols + and x represent the torque on the inner and 
outer surfaces of a virtual cylinder in the fluid, respectively. 
Solid lines are obtained by Eq. (|14[) . The torques for the 
inner and outer surfaces have opposite signs; for clarity, only 
the absolute values are shown. Error bars are smaller than 
the size of the symbols. 



MPC-AT— a shows good agreement with a deviation of 
only about 1% for n = 10 (see the inset of Fig. [J). 



D. Parameters 

We simulate two-dimensional flows. The simulation 
data are displayed with the units of length a, time r = 
a \J mo/kB®, and viscosity 770 = V^feQ/a- We use n = 
10 and At/r = 0.05 or 0.1. Since our aim is to clarify the 
difference of viscous stresses between +a and —a fluids, 
we use small angular velocities fir = 0.004 to 0.01 for 
circular Couette flows to keep the density constant and 
a low Reynolds number Re — pD 2 Q/r] 1, where D = 
10a is the diameter of the smaller cylinder. To obtain 
the hydrodynamics of liquids, we use a small Knudsen 
number Kn = l\/D = 0.01, where l\ = At^Jk^Q/rno is 
the mean free path of fluid particles. The error bars are 
estimated from three or ten independent runs. 

MPC-AT is more time consuming than MPC-SR due 
to the heavier use of random numbers (dN c Gaussian- 
distributed instead of d — 1 uniformly-distributed ran- 
dom numbers). On the other hand, taking angular- 
momentum conservation into account only slightly in- 
creases the required CPU time in two dimensional simu- 
lations. 



IV. CIRCULAR COUETTE FLOW 

We consider Couette flow, since it is a well analyzed, 
simple system. Let R\ and R 2 be the radii of two coaxial 
cylinders rotating with the angular frequencies Qi and 



Q2 respectively, where the indices 1 and 2 refer to the 
inner and outer cylinders, respectively. We assume both 
cylinders to be of infinite length and their angular veloc- 
ities to be sufficiently low, such that no Taylor-Couette 
instabilities occur, and the problem can be considered in 
two dimensions. For symmetry reasons, the radial veloc- 
ity component vanishes and the Navier-Stokes equation 
yields the azimuthal velocity [38| 



vg(r) = Ar + B/r 



where 
A = 



£l 2 R 2 
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R\ 
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The torques acting on the cylinders in an +a fluid, which 
conserves angular momentum are [38| 



Tx = -T 2 = 



R{ 



(13) 



The torque T in a fluid at radius R\ < r < R 2 is calcu- 
lated from the momentum transfer across a virtual cylin- 
der of radius R = r, and is equal to T\ or —T\ on the 
inner or outer surface of the virtual cylinder, respectively. 
Thus, the torque on the inner cylinder propagates to the 
outer cylinder via the fluid with a constant value because 
of angular-momentum conservation. However, in an —a 
fluid, the negligence of the angular-momentum conserva- 
tion generates an additional torque. 



A. Uniform Angular Velocity 

First, we consider the simplest case, where the whole 
fluid rotates with constant angular velocity fio- This is 
done with R\ = and fl 2 = or both cylinders rotate 
with the same angular velocity Qi = Q 2 = fio- Here, 
no torque is expected to be acting on the cylinders in 
+a fluids, as this corresponds to the rotation of a rigid 
body. The MPC-AT+a simulations yield the physically 
correct result, T = 0, at any r. However, in the MPC- 
AT— a and MPC-SR simulations, we do observe positive 
or negative torques on the confining inner (R = Ri) or 
outer (R = R 2 ) cylinder, respectively. In the following, 
we consider the torques on the inner and outer surfaces 
of a virtual cylinder of radius R in the fluid, which shows 
the torque generation in the —a fluid. In the MPC sim- 
ulations, we calculate the torques on the inner and outer 
surface of this virtual cylinder by measuring the change 
of the angular momentum per time step in cells crossing 
the virtual cylinder at R = r. The results are shown 
in Fig. [21 This torque is explained by the stress term 
of the asymmetric viscosity f\ in Eq. ([3]). The torque is 
the tangential stress 2t)Q,q multiplied by the circumfer- 
ence length 27r i? and the radius R, i. e. \T\ = AnfjVlQR 2 . 
The torque T m — (|Tj n | + |T out |)/2 averaged on inner and 
outer surfaces agrees with this prediction. 
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However, inner and outer surfaces of the cylinder re- 
ceive slightly smaller and larger torques than T m . This 
mismatch is qualitatively explained as follows. The total 
transferred momentum of particles inside the cylinder is 
equal to that of outer particles with the opposite sign, 
since the translational momentum is conserved. Thus, 
the torque of inner particles is smaller than the outer 
one, since the average distance from the cylinder axis 
of inner particles is smaller. In order to calculate this 
finite-cell-size effect quantitatively, we consider the trans- 
fer of momentum crossing a cylinder of radius R in the 
fluid. It is derived in analogy to the momentum crossing 
a plane in calculations of the viscosity [H, [3, [TH, HU , and 
the details are described in the Appendix. The resulting 
torques T; n and T ou t of a virtual cylinder of radius R that 
are exerted on the inner and outer surfaces, respectively, 
are found to be 



T, ,,(/>') t-l;7,/!W 2 ( 1 : ^j. 



(14) 



Thus, the first-order correction term is +3a/4i?. The 
same correction term can also be derived for —a versions 
of the other MPC methods. This correction term well 
describes the torque difference between inner and outer 
surfaces (see Fig. [5]). 




B. Angular Velocity Gradient 

Next, we consider the flow with angular velocity gra- 
dient induced by (Sli , r2 2 ) = (0,Slo) or (Sl , 0). Both 
+a and —a fluids yield the velocity field described by 
Eqs. CP and (HU) [see Fig. G^a)]. The torque T in the 
+a fluid is constant throughout the fluid, and depends 
only on the relative angular velocity Q± — SI2, not on the 
absolute value of Sli or f^. This results agrees with the 
prediction of Eq. (fT3|) . However, the torque T in the —a 
fluid is not constant and depends on the value of the an- 
gular velocity because of the non-conservation of angular 
momentum. This dependence is well described by Eq. ([3]) 
[see Fig. Hb)]. 

The inset of Fig. [3)Ja) shows the slip velocity Av$ = 
v g sl — Vg h on the boundaries, where Vg h is given by 
Eqs. (fTTj) and (|12p. The velocity v g st is calculated from 
a least-squares fit to Eq. (fTTj) with parameters A and 
B for the range 6 < r/a < 9. The 'u-gra' algorithm 
shows very small slip and the velocity in Fig. [3ja) coin- 
cides with the theoretical values very well. The 'cj-con' 
and 'cen' algorithms show larger slip and ±a fluids show 
similar dependence. 

C. Phase-Separated Binary Fluids 

A boundary of a fluid exists not only on solid objects 
but also between two fluids or on membranes. In order to 
investigate the fluid-fluid boundary in —a fluids, we con- 
sider binary fluids with a fixed geometry of the boundary 



FIG. 3: (Color online) (a) Azimuthal velocity vg and (b) av- 
eraged torque T m = [7j n | + |T ou t|)/2 of circular Couette flow 
for n = 10 and At = 0.1. The inner (Ri = 5a) or outer 
(i?2 = 10a) cylinder rotates with £Iqt = 0.01, the other cylin- 
der is fixed (SI = 0). (a) Symbols represent MPC-AT+a (o, □) 
and MPC-AT-a (+, x). (b) Symbols represent MPC-AT+a 
(•, x) and MPC-AT-a (o, □). Solid lines are obtained by 
(a) Eqs. 113, JlJl) and (b) Eqs. ([I3j), ©. Error bars are 
smaller than the size of the symbols. The inset shows the slip 
velocity Ave (+a, —a) on inner (□, x) and outer (o, A) cylin- 
ders in the outer cylinder rotation for the different boundary 
algorithms introduced in Sec. MI Bl 



surface, which is impenetrable to the fluid particles. The 
inner cylinder of radius R\ of circular Couette flow is re- 
placed by a more viscous fluid, and the outer cylinder 
with radius i?2 rotates with constant velocity Q2 = f2o- 
This is a simplified description of oil and water phase- 
separated due to surface tension, or two liquids sepa- 
rated by a membrane. It is assumed that cylinders ro- 
tate very slowly, and that the flow stress does not change 
the shape of the interface. In MPC- AT, the fluids inside 
(r < R\ = 5a) and outside (R\ < r < R 2 — 10a) have 
high viscosity 771 with mass mi and low viscosity 772 with 
mass mo, respectively (note that r\ cx m). The particles 
of both fluids are scattered elastically at the boundary 
surface at R\ during the streaming step, but the MPC 
collision performed in cells crossing the boundary prop- 
agates the momentum from one fluid to the other. 

In MPC-AT+a, both fluids rotate with S1q indepen- 
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FIG. 4: (Color online) Azimuthal velocity of binary fluids in 
a rotating cylinder with Qo ~ 0.01/r. The viscous fluids with 
particle mass mi and mo are located at r < 5a and 5a < r < 
10a, respectively. Symbols represent the simulation results of 
MPC-AT— a with mi /mo = 2 (+) or mi /mo = 5 (x), and 
MPC-AT+a for mi/mo = 5 (o). Solid lines represent the 
analytical results for MPC-AT— a at mi /mo = 5. Error bars 
are smaller than the size of the symbols. 

dent of their viscosities. However, in MPC-AT— a, the 
inner fluid rotates more slowly for mi > too (see Fig. |4}. 
This is caused by the asymmetric stress term 2f)Q for —a 
fluids where fj ~ r) co \/2. If both fluids rotate at the same 
angular velocity, the inner and outer stresses do not co- 
incide. Thus, the angular velocity of the inner fluid fl\ is 
smaller than the outer one. The inner and outer flows are 
described by Vg(r) = flir and and Eq. (fTTj) . respectively. 
Then, Qi is obtained from the stress balance at r = R±, 
i. e. 2?7if2i = (8/3)772(^0 — ^1) + 2772^1 . This calculation 
well reproduces the numerical results (see Fig. 2]). Thus, 
it is essential to employ an +a version of MPC in simula- 
tions of multi-phase flows of binary fluids with different 
viscosities. 



D. Ideal Star Polymers 

MPC simulations have been used intensively to inves- 
tigate the behavior of macromolecules under flow [1, [H, 
[20 . 12H [22l . I23I [23 | . Here, we consider a two-dimensional 
ideal star polymer with / arms and arm length Lj in a 
MPC fluid, where the central monomer is fixed in the 
center of the enclosing cylinder with Ri — 10a, which 
rotates with constant angular velocity Qq. Consecu- 
tive monomers are connected by the harmonic potential 
W n — §(r„ — iVi+1) 2 but are otherwise not interacting 
with each other. The coupling to the fluid is achieved 
by including the monomers of mass M in the collision 
step [lU . We simulate stars with / = 5 and / = 10 
arms, both with an arm length Lf — 10. We choose 
M = 5to and K — 2k-QT/a 2 for the spring constant, 
i. e. ((r„ — r n+ i) 2 ) / a 2 = 1 in equilibrium. 

We determine the average azimuthal velocity of fluid 




FIG. 5: (Color online) Azimuthal velocity ve of a fluid with a 
star polymer fixed in the center of a rotating cylinder. Circles 
represent the simulation results for MPC-AT+a, whereas the 
results for MPC-AT— a are shown as squares, with full squares 
for / = 10 and open squares for / = 5. The line shows the 
theoretical result for an angular-momentum conserving fluid. 
The inset shows the corresponding radial monomer number 
density distributions for / = 10 (full line) and / = 5 (dashed 
line). 

particles and monomers as a function of r, as shown 
in Fig. El While the MPC-AT+a yields the physically 
correct result, we find a non-uniform angular velocity in 
MPC-AT— a fluid, similar to the case of the binary fluid, 
but without a sharp interface. The star polymer, which is 
located at small radii (see density distribution in the inset 
of Fig. O, rotates more slowly than the cylinder with an 
average angular velocity r2 s tar/^o = 0.63+0.01 for / = 10 
and O s tar/^o = 0.74 ± 0.01 for / = 5. Note that for the 
chosen parameters, the radial monomer density near the 
center is quite large, see the inset of Fig. O The effect 
of a reduced angular velocity is less pronounced for less 
compact stars, i. e. for reduced arm number or decreased 
spring constant. Thus, this artifact can be drastically 
reduced by keeping the local monomer density low, for 
example by taking into account excluded volume interac- 
tions. The —a methods should not be employed for high 
local density of embedded objects. 

V. ECCENTRIC CYLINDERS 

Going one step further, we study a fluid between eccen- 
tric cylinders with radii R± and R2 and fixed axes. The 
outer cylinder is stationary and the inner one is rotating 
about its axis with an angular velocity fii (see Fig. [6]). 
Neglecting inertial forces, Muller [39j derived theoretical 
expressions for the arising torques and forces acting on 
the cylinders, where the latter is predicted to be perpen- 
dicular to the line connecting the two centers. We per- 
form MPC-AT+a simulations with constant torque, and 
measure the resulting angular velocities and the forces 
acting on the inner cylinder as a function of the axis offset 
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d. To avoid any bias, we use the simple 'cen'-boundary 
condition (see Sec. IIII B|) . The results are compared in 
Figs. [7] and [8J with the theoretical predictions of Ref. 
[39(. In general, good agreement is found, although the 
cylinder is rotating up to 7% faster than theoretically ex- 
pected. This can be explained by the finite slip on the 
surface of the cylinder; in order to suppress the slip com- 
pletely, an extrapolation of the velocity field would be 
necessary for the virtual particles, as discussed in Sec. 
IIII Bl Moreover, we also observe a small radial compo- 
nent of the force, shown in the inset of Fig. [51 It tends to 
move the inner cylinder to the center of the outer one, as 
it is expected when inertial effects are taken into account 
401. 




0.4 0.6 
d/(R 2 -Ri) 




FIG. 6: Geometry a fluid (gray) enclosed between eccentric 
cylinders. 



0.007 




0.006 



3-0.005 



0.004 



0.003 



FIG. 7: Angular velocity £l\ as a function of the axis offset 
d for a constant torque T = 75mor _2 a 2 . The squares are 
the results for the MPC-AT+a simulation compared to the 
theoretical result of Ref. [39;] (full line) , the dashed line serves 
as a guide to the eye. The radii of the cylinders are Ri = 10a 
and R2 = 20a. The parameters for the fluid are n — 10 and 
At — 0.05r. Error bars are smaller than the symbol size and 
are therefore omitted. 



FIG. 8: Forces F x , F y on the rotating inner cylinder perpen- 
dicular and parallel to the line connecting the centers of the 
cylinders. The squares are simulation data obtained from the 
MPC-AT+a simulations, compared to the theory of Ref. [39l ] 
(full line). 



In the following we consider one specific geometry in 
more detail, depicted in Fig. [6] with R\ — 10a, R2 — 
20a and an axis offset d = 8a, using the MPC-AT+a 
algorithm. 




FIG. 9: Force distribution on the surface of the counter- 
clockwise rotating inner cylinder (indicated by the dashed 
line) for the geometry of Fig. [6] (axis offset d = 8a) . A con- 



stant torque T = 75mor 



is applied to the inner cylin- 



der, resulting in a counter-clockwise rotation with (fii) = 
0.0035-r" 1 . The length of the arrows is proportional to the 
local force. The direction of the resulting total force is marked 
by the thick arrow acting on the center of the cylinder. 

In Fig. [9[ we present the measured force distribution 
on the surface of the inner cylinder, where the force due 
to the isotropic hydrostatic pressure has been subtracted. 
The force exerted by the fluid on the cylinder is composed 
of two contributions: First of all, the shear stress is hin- 
dering the counter-clockwise rotation. This force tangen- 
tial to the surface is more pronounced in the small-gap 
region (top of Fig. [5] and [5]) than in the large-gap re- 



8 



gion, thus the net force due to viscous stress points in 
the positive x direction, i. e. to the right in Fig. [51 Sec- 
ond, where the fluid is moving into and out of the slit, 
the dynamic pressure gives rise to regions of increased 
pressure on the right side and decreased pressure on the 
left (see Fig.[T0|). This in turn induces a force pointing in 
the negative x direction, counteracting the force due to 
viscous stress and exceeding the latter in strength, hence 
the total force points to the left. Since the MPC fluid is 
compressible, density inhomogeneities emerge, but they 
are sufficiently small that the corresponding local varia- 
tion of the viscosity is negligible. The density distribu- 
tion, which is proportional to the pressure distribution, 
is shown in Fig. [lOl The corresponding stream lines are 
shown in Fig. 1111 resembling very much the theoretical 
results of Refs. [40, |41| . In particular, a back-flow occurs 
in the large-gap region. 




CO 
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FIG. 11: Stream lines for the same parameters as in Fig. [9] 



sity inhomogeneity is slightly less pronounced for MPC- 
AT— a, giving rise to a smaller pressure gradient. This 
is also reflected in the smaller total force acting on the 
rotating inner cylinder in the MPC-AT— a simulation: 
For the MPC-AT-a method we find F x = (-6.16 ± 
0.011)m o ar- 2 compared to F x = (-6.76±0.014)m o ar- 2 
in the MPC-AT+a simulation. 



VI. SUMMARY 



FIG. 10: (Color online) Density distribution of a MPC-AT+a 
fluid enclosed by eccentric cylinders with n = 10 and the same 
parameters as in Fig. [9] 

Next, we study how the results are affected by the lack 
of angular-momentum conservation. Clearly, for a given 
external torque, in the MPC-AT— a method the inner 
cylinder would rotate with an incorrect angular velocity. 
As this effect has already been discussed in the preceding 
section, we fix the velocity on the boundary instead of im- 
posing a constant torque. We choose the results for the 
angular velocity obtained from the MPC-AT+a simula- 
tions as an input parameter for the MPC-AT— a simula- 
tions in order to investigate the influence of the angular- 
momentum conservation on the resulting velocity field. 
For comparability, we chose the parameters in such a 
way that the viscosity r\ — fj + 77 is the same for both 
simulation methods. We find practically identical veloc- 
ity fields. In order to quantify the difference, we calcu- 
late the ratio (( v AT + a - v SR ) 2 ) 1 / 2 /((v AT+a ) 2 ) 1 /2 ~ .03, 
where v AT+a and v SR denote the velocity fields obtained 
by the two different simulation methods, and the average 
is taken over the simulation box. 

Although the density distributions for both simula- 
tion methods are qualitatively very similar, the den- 



We have investigated the relevance of angular- 
momentum conservation in mesoscale hydrodynamics 
simulations. We have focused on MPC methods, but 
similar results are also expected in other hydrodynamic 
methods without angular-momentum conservation, such 
as DPD— a [29]. Focusing on fluids confined between ro- 
tating cylinders, we compare two simulation variants that 
only differ in the conservation of angular momentum. 

In the bulk, both simulation methods show physically 
correct flow behavior. Here, the negligence of angular- 
momentum conservation simply leads to a modified vis- 
cosity. However, we find that in situations where torques 
are acting on surfaces, often quantitative or even qual- 
itative incorrect results arc obtained without angular- 
momentum conservation. In particular, there are non- 
physical torques occurring even in the rigid body rota- 
tion where no torque is expected. This can be well un- 
derstood from basic continuum fluid mechanics as the 
non-conservation of angular momentum gives rise to an 
asymmetric stress tensor. 

The angular-momentum conservation is essential to be 
taken into account in the following cases to avoid non- 
physical torques, (i) The boundary condition on walls is 
given by forces including torques, such as in circular Cou- 
ette flow, (ii) Finite-sized objects with angular degrees of 
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freedom, or densely distributed point-like objects, rotate 
in fluids by the hydrodynamic stress, such as in colloidal 
and polymer suspensions, (iii) Fluids with different vis- 
cosities are in contact. When the boundary conditions 
are given by velocities, —a methods give the correct veloc- 
ity field. For example, MPC-SR reproduces the frequency 
of von Karman vortex shedding observed in experiments 
and other numerical methods (see the Strouhal number 
in Fig. 5 of Ref. [ll|). Thus, the +a version of mesoscale 
hydrodynamics methods have to be employed whenever 
torques play a role in the flow of (complex) fluids. 
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APPENDIX: Calculation of Torque on Cylinder 
Surfaces 

In order to calculate the momentum crossing a cylinder 
of radius R, we use the following simplification: Instead 
of explicitly taking into account all quadratic cells that 
are intersected by R, due to the cylindrical symmetry 
of the problem it is favorable to consider an annular arc 
with radial width a and area a9 c R c — a 2 as an 'adapted' 
collision cell, where 9 C is the angular width and the inner 
or outer radius is R c — a/2 or R c + a/2, respectively. The 
collision step locally equalizes the velocity within the cell 
on average, hence the pre-collisional velocity distribution 
v(r) = Qreg is converted into the average azimuthal ve- 
locity 



R c +a/2 6 c /2 



v c e(Rc) 



dr 



dO 



Ur 2 cos(0) 



R c -a/2 -6 c /2 
( " 2 

VLR C 1 



24i? c 



(15) 



This collision accelerates or decelerates the fluid parti- 
cles inside or outside of a virtual cylinder with radius R 



in fluids, respectively. We now calculate the change of 
angular momentum AL illiOUt (i?, R c ) = L[ n out (R, R c ) - 
Lin,out(R, Rc) caused by this alteration of the velocity 
distribution, where the subscript 'in' or 'out' denotes 
the inner or outer surfaces of a cylinder of radius R and 
iin,out or L[ n out denotes the angular momenta before or 
after the collision step, respectively. For the inner sub- 
annulus, we find 



R 



AL m (R,R c ) = J dr 



2nm(n — l)r 2 



R c -a/2 

and analogously for the outer sub-annulus 



v c e (R c ) ~ v e (r)} 
(16) 



Rc+a/2 



AL out (R, R c 



dr 



2irm(n — l)r 2 



[v c e (R c )-v e (r)}. 



(17) 

In the derivation, (v$ — v^) = (1— l/n)((vj)— Vg(R c )ee) is 
employed. Note that AL[ n (R,R c ) and AL out (R,R c ) are 
not exactly oppositely equal, reflecting the fact that the 
total angular momentum of the considered annuli slightly 
changes. 

To take into account the random grid shift, we subse- 
quently average over all annuli containing R, i. e. R — 
a/2 < R c < R + a/2. Finally, the torque T in . out = 
(AL in . out (R, R c )) / At is given by 



At 



± 



5R 2 
36 



aR 29a 2 
~~8~ 1080 



R" 



— ± — In 

36 V 8 a J 2R-a 



2R + a 



■nm(n - 1)QR 2 
6At 



±1-^ 
AR 



47r?7fii? z ±1 



3 a 
AR 



(18) 
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